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Primordial black holes (PBHs) are black holes which may have formed very early on during the 
radiation dominated era in the early universe. We present here a method by which the large scale 
perturbations in the density of primordial black holes may be used to place tight constraints on non- 
gaussianity if PBHs account for dark matter (DM). The presence of local-type non-gaussianity is 
known to have a significant effect on the abundance of primordial black holes, and modal coupling 
from the observed CMB scale modes can significantly alter the number density of PBHs that 
form within different regions of the universe, which appear as DM isocurvature modes. Using 
the recent Planck constraints on isocurvature perturbations, we show that PBHs are excluded as 
DM candidates for even very small local-type non-gaussianity, |/tvl| ~ 0.001 and remarkably the 
constraint on g^L is almost as strong. Even small non-gaussianity is excluded if DM is composed 
of PBHs. If local non-Gaussianity is ever detected on CMB scales, the constraints on the fraction 
of the universe collapsing into PBHs (which are massive enough to have not yet evaporated) will 
become much tighter. 
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1 Introduction 

Primordial black holes (PBHs) are black holes which theoretical arguments suggest might have 
formed from the direct gravitational collapse of large density perturbations very shortly after the 
end of inflation. PBHs may theoretically form with any mass, although their abundance is typically 
well constrained by observations. Whilst PBHs with mass lower than 10 15 g would have evaporated 
by today (with the possible exception of Planck mass relics), more massive PBHs would still survive, 
and represent a viable dark matter (DM) candidate. 

Many efforts have been made to observe PBHs, and whilst they have not yet been seen, this 
has led to many corresponding constraints on their abundance in different mass ranges pQ. The 
econstraints typically assume that PBHs form at a single mass scale and are stated in terms of the 
mass fraction of the universe going into PBHs at the time of formation, (3. There exists only a 
narrow window in which PBHs of a single mass could make up the entirety of DM, with other scales 
being excluded by observations. It is noted that there has been a recent claim that the tidal capture 
of PBHs by neutron stars could be used to exclude the remaining window (apart from Planck mass 
remnants) [2], but this has been refuted in [3J 0]- The results presented here can also be applied 
if DM is composed of smaller PBHs which have all but evaporated by today leaving Planck mass 
remnants which may make up DM |5|. Whilst this mass range is not explicitly considered, it is 
certainly not ruled out by observations, and the results presented here are almost independent of 
the PBH mass. 
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In order for a significant number of PBHs to form, the power spectrum on small scales needs to be 
significantly larger than observed in the CMB - of order 10 -2 in the case of gaussian perturbations. 
This is possible in many models of inflation, including the running mass model [6], axion inflation 
0, a waterfall transition during hybrid inflation 0 0 HD] , from passive density fluctuations mi, 
or in inflationary models with small field excursions but which are tuned to produce a large tensor- 
to-scalar ratio on large scales m ■ See also mum EE and a summary of various models which 
can produce PBHs is presented in M- Alternatively, the constraint on the formation criteria can 
be relaxed during a phase transition in the early universe, causing PBHs to form preferentially at 
that mass scale m - although such an effect will not be considered here. 

PBHs have traditionally been used to investigate the early universe by placing a constraint on 
the small scale power spectrum from the corresponding constraint on their abundance [IBlESlEo]. 
In this paper, large scale fluctuations in the PBH density caused by local-type non-gaussianity are 
considered. If DM is composed entirely, or partially, of PBHs, these perturbations will be seen as 
isocurvature modes in cold dark matter (CDM) - upon which there are tight constraints from the 
recent Planck data release m- 

The isocurvature perturbations are formed in a highly non-linear manner in this model. PBHs 
form shortly after horizon reentry during radiation domination, with an energy density exponen¬ 
tially sensitive to the amplitude of the power spectrum. Observational constraints imply that at 
most one region in a million collapsed into a PBH so the large scale radiation density is almost 
unaffected, but if PBHs form DM then the amplitude of the DM perturbation is extremely sensitive 
to the modal coupling. Using this mechanism, CDM (with zero pressure) is formed in a universe 
which could have previously have been made up entirely of radiation and hence had no isocurva¬ 
ture perturbation prior to PBH formation. Such an effect is impossible within linear perturbation 
theory [22] . 

In a previous paper, the peak-background split was used to investigate the effect of modal 
coupling on the constraints which can be placed on the small scale power spectrum [23] . In this 
paper we use the same mechanism to investigate the extent to which modal coupling produces CDM 
isocurvature modes and discuss the implications of such an effect. Even if the initial conditions are 
adiabatic, which has been shown to be the case in single-field inflation, if there is modal coupling 
then the conversion of radiation into CDM (by collapse into PBHs) can have different efficiencies in 
different regions of the universe, which introduces isocurvature modes in the CMD after inflation 
has ended. 

Even single-field inflation generates a small value of /nl with magnitude comparable to the 
spectral index [24] - which apparently could therefore rule out single-field inflation as a mechanism 
to create PBH DM. However, it has been argued that this is a result of gauge choice [25] [26], and 
that for our purposes the effective /nl = 0 in single-field inflation. It is therefore assumed in this 
paper that /nl can be arbitrarily close to zero. 

Throughout, we will assume /tvl to be scale invariant whilst the power spectrum becomes 
several orders of magnitude larger at small scales - which is likely to be unrealistic given a specific 
model. However, this is a conservative approach, because if \/nl\ were to become larger at some 
small scale, it would not weaken the constraints derived here, but would be likely to strengthen 
them. Even if the bispectrum was exactly zero when all three modes have sub CMB scales, the 
modal coupling between the CMB and PBH scales would still effect the amplitude of the power 
spectrum on PBH scales and the constraints which we derive would not be significantly weakened. 
In such a case, the perturbations within a region smaller than we can probe on the CMB would be 
Gaussian, but the variance would vary between different patches, in a way completely correlated 
to the long wavelength perturbation. 

Shortly prior to the release of this paper, Tada and Yokoyama [27] released a paper discussing 


2 


a similar effect and the use of PBHs as biased tracers. We confirm their results and extend 
the calculation to account for the non-gaussianity parameter g^L as well as Jnl, the effect of 
intermediate modes (between the CMB- and PBH-scales), and make use of the more recent results 
from the Planck 2015 data release. Because all surviving PBHs necessarily behave as at least 
a subdominant DM component today, we also show how the allowed fraction of PBHs can be 
constrained more tightly than previously realised, under the presence of even small non-Gaussianity. 

The layout of this paper is as follows: in section 2, the calculation of the PBH abundance, 
in both the gaussian and non-gaussian case, is reviewed. In section 3, modal coupling and how 
the peak-background split may be used to investigate its effects on PBH abundance is discussed. 
In section 4, the calculation is applied to the formation of CDM isocurvature modes and place 
constraints on the non-gaussianity parameters in the case of PBH DM, and the calculation is 
extended to include simultaneous /jvl and gNL, intermediate modes, and the case where PBHs 
only make up a portion of the DM. We conclude with a summary of our arguments in section 6. 

2 Calculating the abundance of PBHs 

The abundance of PBHs is normally stated in terms of j3: the energy fraction of the universe 
going into PBHs at the time of formation. The standard calculation used in the literature uses 
a Press-Schechter approach, although it has been shown that, for a gaussian distribution, this 
matches well when the theory of peaks is used. It has been argued that the density contrast, rather 
than the curvature perturbation, should be used - although an approximation using the curvature 
perturbation works very well if care is taken to exclude super-horizon modes from the calculation, 
and this simplifies the calculation greatly. In this section, we will briefly review the calculation, as 
well as the main sources of error, for both gaussian and non-gaussian cases. 

When a perturbation reenters the horizon, if its amplitude exceeds a certain threshold, or 
critical, value, then gravitational forces will overcome pressure forces and the region will collapse 
to form a primordial black hole. There has been extensive research to calculate the threshold value 
[281 EU EDI EH E2 EH, which is typically stated in terms of the density contrast. The critical 
value of the density perturbation is believed to be 5 C « 0.45. However, in this paper the curvature 
perturbation is used, and the corresponding critical value is ( c « 1 - within the range found by [341 . 
and is consistent with using the density contrast [35] . 

The main source of uncertainty in the critical value is due to the unknown shape of primordial 
perturbations - and this is the largest source of error in the calculation of the abundance. However, 
whilst the effect on the calculated value of the abundance is large, the effect of this uncertainty on 
derived parameters is relatively small. For example, an error of 010% in the threshold value results 
in an error of several orders of magnitude in the calculated (3 but only an error of 010% in the 
constraint on the power spectrum [36], [37] . In this paper, because our results depend only on the 
relative abundance of PBHs in different regions of the universe, the conclusions are not sensitive to 
small changes in the threshold value. 

Using a Press-Schechter approach, the mass fraction of the universe going into PBHs at the 
time of formation is given by integrating over the probability density function (PDF), 

OO 

P = J P(()d(. (1) 

C C 
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In the case of a gaussian distribution, the probability density function is 

P(c) = 75 ? exp (^)' <2) 

Where a 2 is the variance of perturbation amplitude at the PBH forming scale. /3 can therefore be 
written in terms of the complimentary error, 


P 



(3) 


Expanding using the large-x limit of erfc(x), gives 


P 



(4) 


This is valid only if the distribution is gaussian, and because PBHs form in the extreme positive 
tail of the PDF, their abundance is very sensitive to any non-gaussianity, which we discuss below. 


2.1 Calculating the abundance of PBHs in the presence of non-gaussianity 


In the local model of non-gaussianity, the curvature perturbation is given by 

C = Cg + ^J NL (Cg — °' 2 ) + t^9nlCg + •■• = h (Cg) j (5) 

where a 2 is the variance of the gaussian variable Cg> and is subtracted to ensure the expectation 
value of C is zero. 

The calculation of the abundance of PBHs is most easily performed by calculating the values 
of Cg which correspond the critical value, ( c , and integrating over the corresponding regions of the 
gaussian PDF of Cg - the reader is directed to mm for a full derivation. For example, let us 
consider the case where g^L and higher order terms are zero: 


C = Cg + it/jvl (Cg - (j2 ) = h (Cg) • 


h x (Cc) therefore has two solutions, given by 


Kl = h ±(C c) = 


-5± ^25 + OOCc/jvl + 36 Qf 2 NL a 2 


6/atl 
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Furthermore, if we make the assumption that /ml is small, /jvx <C 1, which we will show is justified 
in the case that DM is composed of PBHs (and is further verified by the findings of (27] ), the above 
expressions can be simplified further. In the expression of (3 for positive and negative / nl , the 
first term inside the brackets dominates, and f3 can be written in terms of one complimentary error 
function, 


P = 



d(G 


= erfc 



( 10 ) 


2a 2 


v(h c l ) 2 


exp 


2<J 2 J ' 


Deriving an analytic expression as shown here is not a necessary step, but it is a useful approxima¬ 
tion, and we will later use this result to derive an analytic expression for bias factor and amplitude 
of isocurvature modes in the PBH density. 

Although it is not shown here, the same calculation can be performed for the local model of non- 
gaussianity containing gNL ~ the interested reader is again directed to [36] [37] for a full discussion 
of the calculation. In the case where only a cubic and linear term are considered 

( = (g +—gNL(o = H(g), ( 11 ) 


then /i -1 (C c ) has up to three possible solutions, depending on the value of g^L and 
assuming that g^L is small, g^L <C 1, which again will be shown later, the expression 
by one erfc function as in equation (1101) . with a different expression for h _1 (C c ). To 
9nl 



Cc 


%ldNL 

25 


Cc- However, 
is dominated 
first order in 

( 12 ) 


3 Modal coupling and the peak-background split 


It has previously been shown [35] that curvature perturbation modes which are a long way outside 
the horizon at the time of PBH formation have little effect on whether a PBH forms. This is due to 
the suppression of large scale density modes by a factor k 2 relative to the curvature perturbation. 
In radiation domination: 


6(t, k) 


2(1 + to ) 

5 + 3 w 




(13) 


where u> = 1/3 is the equation of state, and ( aH ) _1 is the horizon scale at the time of PBH 
formation. However, long wavelength modes can have an indirect effect on the abundance of PBHs, 
(3, due to modal coupling from non-gaussianity. A long wavelength mode can affect both the 
amplitude and distribution of the small scale perturbations which may form PBHs. In figure [1] we 
show how the coupling of long- and short-wavelength modes can affect the number of PBHs forming 
in different regions of the universe. At the peak of the long wavelength mode, the amplitude of the 
small scale mode is increased, forming more PBHs, whilst the opposite occurs at the trough. 
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Figure 1: The top plot shows an example of a universe containing only two modes. As an example 
of modal coupling, the amplitude of the short wavelength mode is a function of the long wavelength 
mode - the amplitude of the short-wavelngth mode is larger at the peak of the long-wavelength 
mode. At the time when short-wavelength mode enters the horizon, and PBHs at that scale form, 
the long-wavelength mode is not yet visible and will not affect whether a PBH forms or not. The 
bottom plot shows the same universe, but with the long wavelength mode subtracted, enabling £ 
to be used a formation criterion for PBHs. The dashed red line shows the formation threshold 
for PBHs - regions where the curvature perturbation is greater than the formation threshold will 
collapse to form a PBH. The black circles represent areas which will collapse to form a PBH. It can 
be seen that a relatively small change in the amplitude of the small scale mode can have a large 
impact on the number of PBHs forming in a region. 


How modal coupling can affect the constraints on the power spectrum at small scales from PBHs 
has been investigated [23], although it was assumed that all the modes involved were sub-CMB and 
potentially had a large amplitude. In this paper, we will go beyond previous work and study the 
case where the large-scale modes are observable in the CMB and hence very small. Despite the their 
small amplitude, we show that these perturbations have a remarkably large effect on observations. 
In this section, we will briefly review the calculation using the peak-background split to investigate 
modal coupling due to the local non-gaussianity parameters /nl and gNL, and in the following 
section, apply this to the abundance of PBHs and the creation of isocurvature modes. 

3.1 Quadratic non-gaussianity, 

We will take the model of local non-gaussianity, in terms of the curvature perturbation £, to be 
described by 

C = Cg + jr/iVL (Cg “ °' 2 ) = HCg), (14) 

where (g is a gaussian variable. It is necessary to subtract a 2 = ((q) so that the background 
(average) value of £ remains zero. We will now employ the peak-background split, and write the 
gaussian component as the sum of a long-(background) and short-(peak) wavelength component, 

C G = Cl + Cs ■ (15) 
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Equation (1141) then becomes: 

C = (0 + Cs ) + nl ((Ci + Cs) 2 ~ ((0 + Cs) 2 }) ■ (16) 

However, terms which depend only on the long-wavelength mode do not affect PBH formation, and 
should not be considered when determining the abundance of PBHs. We therefore subtract those 
terms, leaving: 

C = ^1 + jr InlCi ^ Cs + jr (Cs ~ o*) ■ ( 17 ) 

We can now rewrite the expression in terms of new variables, Cg, a and /nl, and calculate the 
abundance of PBHs f3 as described in section 2, as a function of the long wavelength mode, Q. 


C G 


1 + jr InlCi ) Cs, 


a = 

Inl 


1 + jJnlCi ) &s, 


= ( 1 + jJnlCi 


-2 


Inl- 


(18) 


Equation (1171) can then be written in a form analogous to equation (I14|) . 

C = Cg + yJnl (Cg ~ d 2 ) = KCg)- (19) 

Therefore, both the amplitude and distribution of the small-scale perturbations are affected. In 
order to calculate the abundance of PBHs, the variables in equation (|18D can then be inserted into 
equation (fTOl) . 


3.2 Cubic non-gaussianity, g^L 


Here, we will follow the same steps as for /nl, to show how the presence of a cubic term causes 
modal coupling. For this section, we will assume /atl = 0, and C to be given by 

C = Cg + t^9nlCg- ( 20 ) 

Again, using the peak-background split, one obtains: 

C = ^1 + 9NLCl ^ Cs + Cs + Cs + &(Cl), (21) 


where again, the terms dependant only on Q are neglected because they don’t have a significant 
effect on PBH formation. The above expression can then be rewritten in terms of new variables 
Cg, d, f NL and g NL , given by 


Cg = ( 1 + t^9nlCi ^ Cs, 
d = ( 1 + ^SnlCi^J °s, 

Inl = ( ~9nlCiJ f 1 + 7^9nlC? 


-2 


( 22 ) 


. , 27 a2 

9NL = 9NL ( 1 + 25 9nlCi 
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Figure 2: An example of a power spectrum containing a narrow peak. N represents number of e- 
folds, with smaller scales represented by larger N. The power spectrum is small on most scales with 
a spectral index of ns = 0.96, compatible with observations of the cosmic microwave background 
(CMB) and large scale structure (LSS). The narrow peak in the power spectrum corresponds to 
the scale at which will PBHs form. 


Equation (12011 can then be rewritten as 

C = Cg + g/iVL (Cg — + T^SNlCg- (23) 

An expression for the abundance of PBHs in a given region of the universe, /?, can then be derived 
as shown in section 2. 

In this section, it has been shown that long wavelength modes can affect the amplitude of local 
small scale perturbations and the non-gaussianity parameters, and in the next section the effect of 
this on the abundance of PBHs within a given region will be discussed. 


4 The isocurvature modes of PBH DM on CMB scales 


The abundance of PBHs in a region of the universe can be affected significantly by large-scale 
curvature perturbation modes in different regions of the universe. If PBHs make up DM, then 
these differences in the abundance of PBHs will appear as fluctuations in the density of DM. In the 
presence of loacl-type non-gaussianity, the fluctuations in the DM can be significantly greater than 
the curvature perturbations responsible for producing them - and tight constraints can therefore 
be placed on the non-gaussianity parameters if this is the case from the isocurvature constraints 
from Planck. 

We will define the difference in the abundance of PBHs at the time of formation, 5p, as 


8f> = 


tzl 
p ’ 


(24) 


where j3 and j3 are the perturbed and background values of the PBH abundance at the time of 
formation respectively. If the large-scale curvature perturbation ( is small, it can be related to Sp 
by a constant factor b (referred to the scale dependant bias in |27j). 


bp = Kl 


( 25 ) 
















where b is a function of the non-gaussianity parameters, the variance of the small-scale perturbations 
and the critical value for PBH formation £ c . The factor b therefore parameterises the bias of PBHs 
to form in the presence of large-scale curvature perturbations. 

In this section, we will consider the case where the power spectrum is very small on all scales, 
except for a narrow region where there is a sharp spike - which is responsible for the production 
of PBHs of a mass corresponding to this scaled- An example of such a power spectrum is given in 
figure El We therefore ignore in this section the presence of perturbations of intermediate scales, 
but extend the calculation in the following section to account for when there is a broad peak in the 
power spectrum. 

The abundance of PBHs at a later time on a comoving slicing will be affected by difference in 
their density at the time of formation, as well as by the difference in expansion since the time of 
formation - in denser regions of the universe, inflation ends and PBHs form slightly later, so even 
if the PBH density is constant at the time of formation, the density will not be constant. To first 
order in the density of PBHs can be expressed as 

Qpbh = (1 + b( + 3C) &PBHi (26) 


where the 3£ term is simply the adiabatic mode expected from the expansion of the universe, and 
PIpbh is the background density of PBHs. The b(, term therefore represents a deviation from the 
expected amplitude of the mode if it was purely adiabatic - it is an isocurvature mode, which will 
either either be fully correlated, or fully anti-correlated depending on the sign of f nl■ If PBHs 
make up a significant fraction of the DM content of the universe, the constraints on isocurvature 
modes from Planck can then be used to constrain b - and therefore constrain the non-gaussianity 
parameter^). For simplicity in this paper, except section 5.2, we will assume that DM is entirely 
composed of PBHs, and calculate corresponding constraints on the non-gaussianity parameters / a tl 
and gNL■ On CMB scales, the constraints from Planck on isocurvature modes can be used |2T| 


100 / 3^0 


0.13 , fully correlated 
0.08 , fully anti-correlated, 


(27) 


where 


ft. 


Piso 

Piso + Pq 


(28) 


The fully correlated modes correspond to positive b, whilst fully anti-correlated corresponds to 
negative b (and positive/negative /a tl and g^L respectively). The isocurvature power spectrum is 
related to the curvature perturbation power spectrum as 


Piso — b~P^, 


(29) 


and we therefore obtain constraints on b as 


- 0.028 < b < 0.036. (30) 

This result will now be used to derive a result on the non-gaussianity parameters. 

1 The mass of a PBH is roughly equal to the horizon mass at the time of formation. See [35j for further discussion. 
2 Note that the reverse is also true - for a given value of the non-gaussianity parameters, an upper limit can be 
placed on the amount of DM which is made of PBHs 
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4.1 Isocurvature modes from f^L 

In section 2, an expression for the abundance of PBHs at the time of formation /3, was derived 
in terms of the non-gaussianity parameter /jvl, the variance of the gaussian component!! <x 2 , and 
the critical value for collapse Cc - equation (fTOl) . with h given by equation (0. However, this 
calculation assumes there is no coupling to large scale modes (and is equivalent to the background 
value, /?, if large-scale perturbations are small - as is the case here). In section 3 it was shown how 
to account for the presence of a large scale modes - namely, by using the transformed variables /ml 
and a instead, given by equation (fT 8 l) - which calculates the perturbed abundance f3. 

By combing equations 0, (0, Q, (fT 8 l) and it is possible to derive an expression for 8p 
in terms of /nl, &s (where the s subscript has been adopted to denote the small PBH scale), and 
the critical value ( c . Expanding the expression to first order in £ gives the result 

25 + ZOUnl + 36/2 TL a 2 st ~ 5^25 + QOUnl + 36 f* L a* 

< 5/3 =- 7 V -C, (32) 

3 f NL a^ 25 + 60C Jnl + 36 p NL ^ s 

and therefore b is given by 

25 + 30C c/nl + 36/ 2 L a 2 st - 5^25 + GOUnl + 36f^ L a 2 s 

b= - V -, (33) 

3f NL a 2 ^25 + GOUnl + 36 

or to first order in f^L 

b= K 1+ 1) ,kl ' (34) 

As expected, a positive /nl, which boosts the power spectrum on small scales in areas of higher 
density, produces a positive bias, and fully correlated isocurvature modes in PBH DA/@. Negative 
/a tl has the opposite effect, and produces fully anti-correlated isocurvature modes. 

In order to investigate the constraints on the non-gaussianity parameters, it is necessary to 
estimate values for the other parameters involved, and how these would affect the constraints. The 
variance of the small scale perturbations and the critical value. 

• First, Cc is considered: there is significant error in the exact value of the threshold value, due 
to uncertainty in the shape of the primordial perturbation which collapse to form PBHs. Most 
recent simulations have calculated the critical value in terms of the density contrast, finding 
5 C ~ 0.4. This is consistent with the calculation here if the critical value of the curvature 
perturbation is related by a factor |, meaning £ c « 1 , which is consistent with the range of 
values found in figure [34] . Figure [3] shows how the factor b depends on the critical value for 
different values of / nl■ 

3 <j is related to the power spectrum as follows [3HJ 

= <j 2 + Qj (4 fh L + 6pjvi) a 3 4 * In (kL) + 0 j (27 g% L ) a 6 ln(fcL) 2 , (31) 

where the higher order terms from g^L have also been included, and ln(fcL) is a factor of around unity. Note that, 
since the non-gaussianity parameters are found to be very small, the higher order terms will not have a significant 
impact, and to a good approximation V f = a 2 . 

4 The second expression for b corresponds to equation (14) in [27] . The more complicated expression, equation (13311 . 

is because a gaussian distribution on small scales has not been assumed. The differences between the 2 calculations 

are discussed in Appendix B. 
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Figure 3: The plots above show the effects of a different threshold Cc 011 the PBH bias b arising from 
an /tvl term. A larger value of Cc suggests a larger bias factor. The left plot shows the effect for 
negative /nl and the right plot for positive $nl- The dotted black lines represent the constrains on 
b from the constraints on isocurvature modes from Planck. \/nl\ = 8 x 10~ 4 is typically excluded 
whilst |/jvl| = 4 x 10 -4 is typically allowed. To generate these plots the value a = 0.15 has been 
used. 


• To calculate a s , it is necessary to first calculate the value of (3 for which PBHs are otherwise 
unconstrained by observations and could be DM. The range of mass scales in which PBHs 
can form a significant fraction of DM is roughly 10 17 g < Mpbh < 10 24 g [lj. The constraint 
on f3 from the abundance of DM in this range are given by [193 


(3 < 2 x 1(T 19 


/ Mpbh \ 1 ^ 
V/m5 x 10 14 g/ 


(35) 


where /m is the fraction of the horizon mass which ends up inside the PBlfi and Mpbh is the 
mass of the PBH. Assuming DM to be made up entirely of PBHs of a single mass scale within 
this range, f3 can therefore range from (3 < 10 -16 to j3 < 10 -11 . Assuming the most optimistic 
and pessimistic values for (3 and Cc, &s is calculated to lie in the range 0.1 < a s < 0.2 for close 
to gaussian perturbations [36]. Figure [1] displays how b changes with a s . 


Smaller values of the variance of the small scale perturbations, <t 2 , would lead to tighter con¬ 
straints on / nl , whilst a smaller critical value Cc leads to tighter constraints on /nl- Because a 
larger value of Cc implies a larger value of a s , these effects virtually cancel out - and the results 
presented below are therefore not sensitive to uncertainty in Cc- 

Assuming PBH form at a single mass scale, the weakest constraint on Jhl comes from consid¬ 
ering the mass of the largest PBHs which could make up DM, which is taken to be Mpbh = 10 25 g, 
for which f3 ~ 10~ 14 . If DM is made entirely of PBHs, the constraints on f^L are therefore 


-4 x 1(T 4 < f N L < 5 x 1(T 4 . 


(36) 


The results are not significantly different for PBHs of different mass. For example, for Mpbh = 
10 2 °g the constraints on f^L are 


- 3 x 10" 4 < f N L < 4 x 10 


-4 


(37) 


3 /m is a factor of order unity, which is neglected as it has very little effect on the calculated value of a s . 
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Figure 4: The effects of a different a on the PBH bias b arising from an /nl term are investigated. 
A larger value of a suggests a smaller bias factor. The left plot shows the effect for negative /nl 
and the right plot for positive /nl- The dotted black lines represent the constrains on b from 
the constraints on isocurvature modes from Planck. \/nl\ = 7 x 10” 4 is typically excluded whilst 
\/nl\ = 4 x 10 -4 is typically allowed. To generate these plots the value ( c = 1 has been used. 


4.2 Isocurvature modes from g^L 


In addition to /nl, it is interesting to consider isocurvature modes arising from gNL and place 
constraints, or whether the effects of modal coupling from gNL could cancel the effects from /nl- 
The effect of higher order terms are beyond the scope of this paper. 

The same derivation can be followed as that for f nl , leading to an expression for b to first order 


in gNL 


b = 


27 (al 


- C 2 ) K 2 + C 2 ) 

25<7|Cc 


9nl- 


(38) 


Again, as expected, positive gNL corresponds to fully correlated isocurvature modes, and negative 
9nl corresponds to fully anti-correlated isocurvature modes. The PBH bias factor b is again a 
function of the non-gaussianity parameter gNL, the variance of the small scale perturbations ay, 
and the formation threshold Cc- The dependance of b on Cc and a s is shown in figures [5] and [6] 
respectively. 

We see again that smaller values of a s would lead to tighter constraints on gNL, whilst a smaller 
C c leads to tighter constraints on gNL- However, unlike the case with f nl, the constraint which 
can be placed on gNL depends on the value of Cc, although only by a factor of 0(10%). The results 
presented below are the weakest constraints, corresponding to a low formation threshold, for PBHs 
of mass 10 25 g 

- 6 x 10" 4 < g N L < 7 x 10~ 4 . (39) 


Notice that these constraints are very comparable to those on /nl, see (l37l) . The / nl term has 
an effect of 0(1O~ 5 ) on the small-scale power spectrum, whilst the gNL term only has an effect of 
0(1O~ 10 ), and therefore, naively, the constraints on gNL would be expected to be roughly 5 orders 
of magnitude weaker than / nl■ However, a. gNL term also has an effect on the small scale / nl, 
as seen in equation (122(1 . of 0(1O~ 5 ), and because the abundance of PBHs is extremely sensitive 
to non-gaussianity, this causes significant isocurvature modes in the PBH DM. In the case where 
Ci = 10~ 5 and gNL = 10 3 , then /nl ~ 10 -8 . Such a small /nl nonetheless creates a perturbation 
in the PBH density of 0(1C0 6 ), which represents an isocurvature mode of around 10% of C - which is 
excluded by Planck. Because the abundance of PBHs (3 is sensitive to higher order non-gaussianity 
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Figure 5: The plots above show the effects of a different threshold £ c on the PBH bias b arising from 
a 9nl term. A larger value of ( c suggests a larger bias factor. As the expression for b, equation 
(USD, is anti-symmetric under a change of sign of g^L, the results for negative and positive gjvL are 
shown on one plot - but with different constraints on the amplitude of |6|, represented by the dotted 
black lines. \9nlI = 8 x 10~ 4 is typically excluded whilst \9nl\ = 4 x 10~ 4 is typically allowed. To 
generate these plots the value a = 0.15 has been used. 


parameters WL isocurvature modes are expected to rule out significant non-gaussianity at higher 
orders as well - although a quantitative calculation is beyond the scope of this paper. Higher order 
non-gaussianity parameters are considered briefly in section 5.4. 

5 Further consideration of constraints from isocurvature modes 

In section 4, constraints were placed separately on f^L and g^L separately, assuming that DM was 
entirely composed of primordial black holes. In this section, the calculation is extended to account 
for more general models. 

5.1 Isocurvature modes from f^L an d Qnl 

The presence of non-zero non-gaussianity parameters has been shown to create significant isocur¬ 
vature modes, which has led to very tight constraints on these parameters under the assumption 
that DM is composed entirely of PBHs. The calculation is now extended to account for non-zero 
tl and gNL simultaneously - for example, it is possible that the effect of a large positive f^L and 
large negative g^L can cancel out, leaving a very small isocurvature mode. 

Because the non-gaussianity parameters may now become quite large, the full numeric calcu¬ 
lation for the PBH abundance is used to derive a value for the PBH bias b, for example by using 
equations © or Q rather than the much simpler equation (flOl) . 

Figure [7] shows the values of g^L that are permitted for different values of /nl for PBHs of mass 
Mpbh = 10 25 g. Whilst large values of f^L and g^L are allowed, there needs to be significant fine 
tuning to ensure that the resultant isocurvature modes are not excluded by the Planck results - g^L 
needs to have the correct value to 0(0.1%). We note that there is some uncertainty in the value of 
gNL required for a given f^L due to the uncertainty in the formation threshold C, c - although this 
does not affect the conclusion that large non-gaussianity parameters are not allowed unless very 
very finely tuned. This conclusion is expected to remain true for higher-order terms [37j . 
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Figure 6: This plot shows the effects of a different a on the PBH bias b arising from a g^L term. 
A smaller value of a suggests a larger bias factor. As the expression for b , equation (1551) . is anti¬ 
symmetric under a change of sign of qnl, the results for negative and positive gpu, are shown on one 
plot - but with different constraints on the amplitude of |6|, represented by the dotted black lines. 
\9nl\ = 7 x 10” 4 is typically excluded whilst I^wlI = 4 x 10 -4 is typically allowed. To generate 
these plots the value Cc = 1 has been used. This range of a is used because it is approximately 
the range of values required to generate the correct number of PBHs to form DM (assuming that 
perturbations are close to gaussian). 


5.2 Fractional primordial black hole DM 

So far, it has been assumed that DM is made entirely of PBHs. The calculation is now extended to 
account for the fact that PBHs may only make up a small fraction of DM, and this is parameterised 
by fpBH i the ratio of PBH density to DM density. 


rpBH = 


Qpbh 

&DM 


(40) 


In this case, the density of DM is described by 


^dm = (1 + rpBHb( + 3() f 1dm, 


(41) 


and the relative amplitude of the isocurvature modes is now given by rppnb . Therefore, from the 
Planck constraints on isocurvature modes instead give constraints on the factor rppHb, 

- 0.028 < r PBH b < 0.036. (42) 


The constraints which can be placed on the non-gaussianity parameters therefore depend upon 
the PBH DM fraction, rppH • Figure [8] shows the allowed values of Jnl, 9nl and rppH if the PBH 
mass is Mpbh = 10 25 g. 

• Large rppH'- if PBHs make up a large fraction of DM then very tight constraints can be 
placed on the non-gaussianity parameters, /jvl, gNL < 0 (1O~ 2 ). 

• Small rppH'- if PBHs make up a small fraction of DM, rppH < 0.1, then the constraints 
° n Inl and gpn weaken significantly. However, the non-gaussianity parameters only become 
larger than 1 if rppH < 0(1O~ 3 ). In the case where vpbh is very small, the non-gaussianity 
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Figure 7: The constraints on simultaneous f^L and g^L are displayed. The right plot simply 
displays the central region of the plot on the left. The solid lines represent an upper limit from 
fully correlated isocurvature modes, whilst the dotted lines represent a lower limit from fully anti¬ 
correlated isocurvature modes. There is some uncertainty in the value of g^L given a value of /ml 
due to uncertainty in the critical value ( c - the blue lines are obtained using ( c = 0.8, and the red 
lines are obtained using = 1.2. It can nonetheless be seen that large /ml or g^L are excluded 
unless very finely tuned. The shaded regions between the lines can be considered as 2 a contour 
plots from the Planck constraints. 


parameters can become large and it is crucial to account for the effect of a non-gaussian 
distribution on the PBH forming scale, as done in this paper - as seen by the strong asymmetry 
for positive and negative f ml■ 

As tpbh becomes very small, /ml can become large and positive, but is still strongly restricted 
to not be large and negative. This is partly due to the fine tuning of the small scale power spectrum 
necessary to produce a small but not too large number of PBHs when /tvl is negative - even a very 
small amount of modal coupling can mean that this fine tuning is disrupted in different regions 
of the universe, causing large amounts of variation in the number density of PBHs forming. This 
effect is not seen unless the non-gaussian distribution on small scales is accounted for. For gNL, the 
constraints do not depend much on the sign of gNL, and the small difference is due almost entirely 
to the difference in constraints from Planck on fully, or fully anti-, correlated modes. 

5.3 Intermediate modes 

The intermediate scales in between the large scales visible in the CMB and the small scale at which 
PBHs form have so far been ignored. This is a valid approximation if the power spectrum is small 
at all scales except for a narrow peak at the PBH forming scale, as in figure (2j However, this may 
not be the case if, for example, the power spectrum has a broad peak, as seen in figure[9l or becomes 
blue at small scales. In this case, the abundance of PBHs, as well as the amplitude of isocurvature 
modes, can be significantly affected by the presence of perturbations on these intermediate modes. 

If the power spectrum of the intermediate modes is not small, they will have a significant effect 
on the number of PBHs that form, as well as the isocurvature modes visible in PBH DM. This will 
be investigated in a similar to the peak-background split, and the curvature perturbation is split 
into short, intermediate, and long components: 

Cg = Cs + Ci + Cl- (43) 


15 






fPBH 
1.000 r- 

0.500 - 

0.100 - 
0.050 - 

0.010 - 
0.005 - 


>'PBH 




Figure 8: In the case where PBHs only make up a small fraction of the DM content of the universe, 
parameterised by rppH , the constraints on f^L and g^L can become significantly weaker. This is 
due to the fact a large isocurvature mode in the PBH density would only translate into a small 
isocurvature mode in the DM density. The plots above show the allowed values of /jvl and g^L for 
different values of rppH- Whilst the plots show the constraints for PBHs of mass Mpbh = 10 25 g, 
the constraints are not very sensitive to the PBH mass. 


The mass fraction of a given region of the universe going into PBHs is then calculated as before, 
as a function of Q and in addition to / nl, 9nl , cr s and ( c , 

P = HC,iXi). (44) 


However, the intermediate modes are too small scale to be observed in the CMB, and should 
therefore be averaged over: 

OO 

HCl)= f (45) 

—OO 

where /3 is the value of /3 in different (intermediate-scale) regions of the universe, and P(Ci) is the 
probability density function of Q, and is given by: 


P(Ci) = 




(46) 


In principle, (£ 2 ), can be obtained by integrating the power spectrum over the relevant range 
of scales. However, since this is unknown and model dependant, it is parameterised here by r, n t, 
the ratio of the variance of intermediate modes (£ 2 ) to the variance of the short modes u 2 


Tint — 



(47) 


The value of (£ 2 ) can become larger than u 2 due to the fact that many scales can contribute to 
(i, but only one scale contributes to ( s . ((f) is calculated by integrating the power spectrum over 
the range of scales considered to be intermediate 




kmax 



k 

""min 


(48) 
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Figure 9: An example of a power spectrum containing a broad peak. In this paper, there are 3 
difference scales: the large ’’background” scales visible in the CMB, the small ’’peak” PBH forming 
scale (the exact scale of which depends on the mass PBH being considered), and the intermediate 
scales between the background and the peak. In such a case, the intermediate modes can have an 
effect on the PBH bias. 


and can become large if the power spectrum is large over a significant range of this integration. 
In contrast, the PBH scale perturbations are only composed of perturbations from one scal^l. 
Therefore, (Cf) can become significantly larger than aj. even though the power spectrum has its 
largest value at the PBH scale. However, it is likely that in such a scenario, PBHs of multiple mass 
scales would be produced, which is discussed later. 

The amplitude of the isocurvature modes therefore depends on the non-gaussianity parameters, 
the small scale power spectrum < 7 ^, the formation threshold £ c , and r in t . A value for the PBH bias 
b is then calculated numerically, figure [TO] displays b dependant on these variables. The effect of 
intermediate modes on the amplitude of isocurvature modes is relatively small for small /tvl or 
9nl unless the variance of the intermediate scales is very large. The constraints on /tvl can be 
weakened by a factor 0(1), although the constraints on g^L are not significantly affected. 

Note that a model where the power spectrum is large over a broad range of scales would likely 
also produce PBHs with a large range of masses, and vice versa. This fact does not affect the 
conclusions presented here, as the production of PBHs at all mass scales would be affected by bias 
in a similar way. We have shown that intermediate modes can significantly affect the PBH bias, 
although which modes are considered to be intermediate depends on the scale at which PBHs are 
forming, and therefore on the mass of PBHs forming. The exact constraints depend on the form 
of the power spectrum, and must therefore be calculated on a model by model basis, which goes 
beyond the scope of this paper - although the constraints will not be weaker than / nl , 9nl < 10~ 3 . 


5.4 Higher Order terms 

Whilst only the constraints on /tv_l and gjyL have been calculated here, very tight constraints on 
higher order non-gaussianity parameters are also expected. In the same way that a g^L term has 
a small but significant effect on /nl, equation (|5j) , higher order terms affect the previous term. 

6 Formally, a 2 S is given by integrating the power spectrum multiplied by a window function. However, provided 
that the spectral index is close to 1, or alternatively there is a peak spanning approximately 1 e-fold at the PBH 
scale, <Tg is approximately equal to the power spectrum at that scale. 
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Figure 10: The effect of intermediate modes on the PBH bias b is displayed for Mpbh = 10 25 g. 
The variance of the intermediate modes is parameterised by ri n t, the ratio of ((f nt ) to a 2 . The 
effect is negligible unless rj n t becomes large, in which case the PBH bias arising from an /nl term 
becomes significantly smaller, although has little effect for gNL- 


Because the mass fraction of the universe forming PBHs is extremely sensitive to non-gaussianity 
parameters at higher orders m, even very small changes to higher order non-gaussianity parameters 
due to modal coupling creates significant creates significant perturbations in the PBH density at 
large scales. As an example, we will consider a 5 th order term in local-type non-gaussianity: 

81 

C = Cg + ^tpvlCg- ( 49 ) 

Utilising the peak-background split gives a A th order term at small scales, Iinl, given by 

b’NL = (50) 

Inserting Q ~ 10~ 5 and inl = 10 -3 gives Knl ~ 10 -8 . The modulation of the Knl by the long 
wavelength mode Cl then generates a perturbation in the density of PBHs forming, 5p ~ 10~ 6 . In 
the picture of PBH DM, this results in a fully-correlated isocurvature mode, with a bias factor of 
b ~ 0.1 - which is excluded by Planck. Because it can be shown that high order terms have an effect 
on the preceding term which is linear in £, tight constraints are expected on such non-gaussianity 
parameters, only weakening slightly as higher order terms are considered. 

6 Summary 

The effect of modal oupling under the presence of non-gaussianity of the local type produces 
significant isocurvature modes in the density of PBHs in the early universe. If PBHs make up a 
significant fraction of DM, the constraints on isocurvature modes in cold DM from Planck can be 
used to constrain the non-gaussianity parameters - in this paper we have considered /nl and gNL 
Using the constraints from Planck on isocurvature modes enables tight constraints to be placed on 
fNL and gNL, 

I/atl|, \9nl\ < O(10“ 3 ), (51) 

unless /nl and gNL have opposite signs and have been extremely finely tuned so that the effect from 
each term cancels. Cases where the constraints could become weaker have also been considered: if 
the power spectrum is large on scales between those visible in the CMB and the PBH forming scale, 
or if DM is only partially composed of PBHs, finding that under these conditions the constraints 
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weaken very slightly (unless PBHs make up a very tiny fraction of DM). Therefore, the detection of 
significant numbers of PBHs would rule out significant local non-gaussianity, and vice versa. Our 
constraints are almost independent of the PBH mass, and can also be applied to Planck mass relics 
which may be left behind from the evaporation of small PBHs. 

The production of isocurvature modes can therefore be used to constrain PBH forming models 
which may otherwise be permitted. For example, we will consider here two models which may be 
ruled out as mechanisms to produce PBH DM: 

• Hybrid inflation: hybrid inflation typically predicts a non-zero /nl, but there is some freedom 
in the exact value. [39] predicts /ml ~ — 1/iV*, where IV* is the number of e-folds between 
horizon exit of some pivot scale and the end of horizon. Inflation is believed to have lasted at 
least 50 — 60 e-folds, which would give /a jl = 0( 10” 2 ) - several orders of magnitude higher 
than allowed by the constraints presented here. [ID] predicts that /a tl can span a range 
of values from ICR 2 to 10 5 - the entire range of which would be ruled out as a method of 
producing PBH DM. 

• The curvaton: the amount of non-gaussianity in the curvaton model depends on the density 
parameter, Q x , of the curvaton, x, at the time it decays into radiation: f^L = —5/4 if D x = 1 
P~I1] . Although higher order local non-gaussianity terms are generated, it is unlikely that these 
will generate small isocurvature perturbations to evade the constraints. 

There are, however, limitations to the calculations carried out in this paper. Notably, we have 
only considered local-type non-gaussianity, and throughout it has been assumed that /a tl and g^L 
are scale invariant. We have also only calculated the dependance of isocurvature modes on f^L 
and gNLi and shown them to a roughly equivalent effect - with g^L having only a marginally 
smaller effect. Higher order terms are therefore also likely to have a similar effect on isocurvature 
modes. We also note that it has recently been observed that sub-horizon perturbations at the time 
of PBH formation have an effect on whether a perturbation will collapse to form a PBH or not 
[12] . The expected amplitude of these sub-horizon modes would be affected by modal coupling - 
and therefore affect the amount of PBHs forming, affecting the isocurvature modes. However, this 
effect is expected to be negligible whilst the non-gaussianity parameters are very small. 
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A Full expression for 5/g from a g^L term 

For completeness, the full expression for 5p arising from a g^L term is included - though this 
expression is still only valid for small gNL- This expression would replace the simpler equation 

m- 
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b = r-503 1 / 3 10 2 / 3 + (-27+ ^300 + 729 sjvlCc 2 ) 7 (45 S^IO^V^ZC 


-5 3 1 / 6 10 1 / 3 100 + 243 5iVL Cc 2 + 100 (-27 yfg^lCc + ^300 + 729 5jvl C, 


1/3 


-54 g NL a 2 (-27V^ZCc + ^300 + 729^^ ' j J (-253 1 / 3 10 2 / 3 ^ (100 + 243 ff7VL Cc 2 ) 


^VgWICc + ^300 + 729^lCc 2 ) ' ~ 250^3 g% 2 L (V^lO^lOO + 243 5 jvl C c 2 - 20 (-27^^10 


+^/300 + 729 5jvl C c 2 ) 7 J + 225 ^ l Cc (^30 3 1/6 10 1/3 - 6yW+243^C?(-27V^ZCc 


+^300 + 729stvlCc 2 ) + 3 5/6 10 2 / 3 (-27V^ZCc + ^300 + 729 5 jvlC 
+243 3 1 / 6 10 1 / 3 ^ L a 2 C c ^30 + 10 1 / 3 (-SlV^ZCc + 3^300 + 729 5 jvl Cc : 

-27 f—4503 1 / 6 Cc 2 (-27^ZCc + ^300 + 729 57V lCc 2 ) ' 

+10 1 /V-\/l00 + 243 sjvlCc 2 ^10 3 1//3 + 10 1 / 3 (-27V^ZCc + ^300 + 729 5 jvlC 
^150 30 1/3 ^y^- + 243C c 2 (-27 ^ l Cc + V 7 ^ 3 ( 10 ° + 243 <WlCc 2 )) ' (-103 1 / 3 

+10 1 / 3 (-27V^ICc + ^300 + 729 57 vlCc 2 ) 7 j J C 


B Comparison with “Primordial black holes as biased tracers” 

In their paper, ’’Primordial black holes as biased tracers ” m, Tada and Yokoyama derive an 
expression for the scale-dependant bias given by 

A6(fc) = 2 f NL MY\k)%. (52) 

This is equivalent to equation fl34l) in this paper. The factor of 3/5 difference is due to a different 
definition of /nl, and the factor A4^ 1 (fc) is a result of their use of the density contrast rather than 
the curvature perturbation. The +1 in the brackets of equation (|34l) is a small correction and can 
be neglected. Therefore, the results for very small f^L in the 2 papers are equivalent. In figure 
ELD the two expressions are compared. For \/nl\ < 0(10 2 ) the two calculations match well, but 
diverge rapidly for larger |/atl|- 

It is therefore necessary to use the full calculation derived in this paper in situations where 
/nl could become larger than 10~ 2 . Whilst such a large value of / nl is generally excluded by the 
constraints on isocurvature modes in the PBH DM scenario, it is relevant where higher order terms 
are considered, or that PBHs form a sub-dominant component of DM. 
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Figure 11: A comparison of the results derived in this paper with those derived in [27] . The solid 
red line denotes the full expression for the PBH bias given by equation (l33lh and the dashed blue 
line represent the scale-dependant bias given by equation (14) in [27] . To make these plots, the 
values a s = 0.1 and ( c = 1 have been used. 
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